#!/bin/bash

export BASE_DATA_PATH=/home/mixagol/data/
export BLAST_DB_PATH="${BASE_DATA_PATH}/1_databases/blast_db/GENOMES_DB"
export SORT_BIN=/home/mixagol/coreutils-8.13/src/sort

export LC_ALL=C
export GENOMES_NUM=`cat "${BASE_DATA_PATH}/3_normed_matrix/map_int_genom.txt" | wc -l`
g++ -O3 6_annotation/intersect_genes.cpp -o 6_annotation/intersect_genes
g++ -O3 6_annotation/intersect_genes_avg.cpp -o 6_annotation/intersect_genes_avg

export ANNOTATE_DIR=${BASE_DATA_PATH}/6_annotation/annot_$(date +%s)
mkdir -p $ANNOTATE_DIR

cat $1 > $ANNOTATE_DIR/gene.faa

cat $ANNOTATE_DIR/gene.faa | 2_raw_matrix/run_blast.sh $BLAST_DB_PATH > $ANNOTATE_DIR/blast_results.txt

for E_VALUE in 5e-02 1e-03 1e-05 1e-11 1e-20 1e-50; do

    export E_DIR=$ANNOTATE_DIR/e_$E_VALUE
    mkdir -p $E_DIR

    cat $ANNOTATE_DIR/blast_results.txt | awk -F'\t' '{if ($11<='$E_VALUE') print $1"\t"$2}' \
        | 3_normed_matrix/apply_maps.py \
            --genomes-map ${BASE_DATA_PATH}/3_normed_matrix/map_int_genom.txt \
        | 3_normed_matrix/create_matrix.py \
        > $E_DIR/vectors.txt
    
    if [ $(cut -f2 $E_DIR/vectors.txt | wc -l) -gt 0 ]; then
        break
    fi
    
    6_annotation/intersect_genes \
        -m ${BASE_DATA_PATH}/3_normed_matrix/matrix_${E_VALUE}_.txt \
        -v $E_DIR/vectors.txt \
        -g $GENOMES_NUM \
        > $E_DIR/joined_raw.txt
    
    join -t$'\t' -1 1 -2 2 \
        <($SORT_BIN -t$'\t' -k1,1 -S100M ${BASE_DATA_PATH}/3_normed_matrix/map_int_gene.txt) \
        <($SORT_BIN -t$'\t' -k2,2 -S100M ${E_DIR}/joined_raw.txt) \
        | cut -f2,4 \
        | $SORT_BIN -t$'\t' -k1,1 -S100M \
        > ${E_DIR}/joined.txt
    
    export PROB_LIMIT=1e-7
    
    join -t$'\t' \
        ${E_DIR}/joined.txt \
        ${BASE_DATA_PATH}/5_pathways/gene_kegg_map.txt \
        | awk -F'\t' '{if ($2<'"$PROB_LIMIT"') print $0}' \
        > ${E_DIR}/joined_kegg_limited.txt
    
    ##
    ## probability<TAB><gene_kegg_id<TAB>function
    ##
    join -t$'\t' \
        ${E_DIR}/joined_kegg_limited.txt \
        <(cat ${BASE_DATA_PATH}/1_databases/genes_info.txt ) \
        | awk -F'\t' '{if (($6 != "__NULL__") && ($6 != "hypothetical protein") && (index($6, "unknown") != 1)) print $2"\t"$3"\t"$6}' \
        | $SORT_BIN -t$'\t' -g -k1,1 -S100M \
        > ${E_DIR}/results_gene_functions.txt
    cat ${E_DIR}/results_gene_functions.txt | head -n100 > ${E_DIR}/results_gene_functions_top100.txt

    ##
    ## probability<TAB><gene_kegg_id<TAB>path_id<TAB>pathname
    ##
    join -t$'\t' -1 2 \
        <(join -t$'\t' -2 2 \
            ${BASE_DATA_PATH}/5_pathways/kegg_pathway_map_full.txt \
            <(cat ${E_DIR}/results_gene_functions.txt | cut -f1,2 | $SORT_BIN -t$'\t' -k2,2 -S100M) \
        | $SORT_BIN -t$'\t' -k2,2 -S100M \
        ) \
        <(cat ${BASE_DATA_PATH}/5_pathways/pathways_full.txt | cut -f1,4) \
        | awk -F'\t' '{print $3"\t"$1"\t"$2"\t"$4}' \
        | $SORT_BIN -t$'\t' -g -k1,1 -S100M \
        > ${E_DIR}/results_pathways.txt        
        
    join -t$'\t'  \
        <(6_annotation/intersect_genes_avg \
            -m ${BASE_DATA_PATH}/5_pathways/pathway_vectors_avg_${E_VALUE}_.txt \
            -v $E_DIR/vectors.txt \
            -g $GENOMES_NUM \
            | cut -f2- | $SORT_BIN -t$'\t' -k1,1 -S100M) \
        <(cat ~/data/5_pathways/pathways_full.txt | cut -f1,4) \
        | awk -F'\t' '{print $2"\t"$1"\t"$3"\t"$4}' \
        | $SORT_BIN -t$'\t' -g -k1,1 -S100M \
        > $E_DIR/results_pathways_avg.txt
    
done

    ##
    ## группируем гены по путям, сортируем группы по числу генов
    ##
    #join -t$'\t' -2 3 \
    #    <(join -t$'\t' -1 2 \
    #        <(cat ${E_DIR}/results_${GENE_ID}_limited_${PROB_LIMIT}_.txt \
    #            | cut -f4 | sort | uniq -c | awk '{print $1"\t"$2}' ) \
    #        <(cat ${BASE_DATA_PATH}/5_pathways/pathways_full.txt | cut -f1,2,4)
    #    ) \
    #    <(cat ${E_DIR}/results_${GENE_ID}_limited_${PROB_LIMIT}_.txt | cut -f2,3,4,5 | $SORT_BIN -t$'\t' -k3,3 -S100M) \
    #    | $SORT_BIN -t$'\t' -k2,2nr -k1,1 \
    #    | 6_annotation/group_genes_from_pathway.py \
    #    > ${E_DIR}/results_${GENE_ID}_pathways.txt


    #cat $E_DIR/results_pathways_avg.txt \
    #    | awk '{if ($1<0.0000001) print}' \
    #    | cut -f4 \
    #    | $SORT_BIN \
    #    | uniq -c | sort -nr -k1,1 \
    #    > $E_DIR/results_pathways_avg_uniq.txt
